Hi again! We hope you are excited to see how the concepts you just learned about work in practice! Let’s dive right in and go through some examples!
If you’re running Mac or Linux, the previous command to install sf may not work first time. These operating systems (OSs) have ‘systems requirements’ that are described in the package’s README. Various OS-specific instructions can be found online, such as the article Installation of R 4.0 on Ubuntu 20.04 on the blog rtask.thinkr.fr.
#install.packages("sf")
#install.packages("spData")
sf - Simple Features for R
spData - Diverse spatial datasets for demonstrating, benchmarking and teaching spatial data analysis. It includes R data of class sf (defined by the package ‘sf’)
library(sf)
## Linking to GEOS 3.8.1, GDAL 3.2.1, PROJ 7.2.1
library(spData)
## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`
library(ggplot2)
library(patchwork)
The object loaded is a sf object containing a world map data from Natural Earth with a few variables from World Bank
# Alternative 1
world = spData::world
# Alternative 2
world = read_sf(system.file("shapes/world.gpkg", package = "spData"))
colnames(world)
## [1] "iso_a2" "name_long" "continent" "region_un" "subregion" "type"
## [7] "area_km2" "pop" "lifeExp" "gdpPercap" "geom"
head(world$name_long, 20)
## [1] "Fiji" "Tanzania"
## [3] "Western Sahara" "Canada"
## [5] "United States" "Kazakhstan"
## [7] "Uzbekistan" "Papua New Guinea"
## [9] "Indonesia" "Argentina"
## [11] "Chile" "Democratic Republic of the Congo"
## [13] "Somalia" "Kenya"
## [15] "Sudan" "Chad"
## [17] "Haiti" "Dominican Republic"
## [19] "Russian Federation" "Bahamas"
methods(class = "sf")
## [1] [ [[<- $<-
## [4] aggregate as.data.frame cbind
## [7] coerce dbDataType dbWriteTable
## [10] filter identify initialize
## [13] merge plot print
## [16] rbind show slotsFromS3
## [19] st_agr st_agr<- st_area
## [22] st_as_s2 st_as_sf st_bbox
## [25] st_boundary st_buffer st_cast
## [28] st_centroid st_collection_extract st_convex_hull
## [31] st_coordinates st_crop st_crs
## [34] st_crs<- st_difference st_filter
## [37] st_geometry st_geometry<- st_inscribed_circle
## [40] st_interpolate_aw st_intersection st_intersects
## [43] st_is_valid st_is st_join
## [46] st_line_merge st_m_range st_make_valid
## [49] st_nearest_points st_node st_normalize
## [52] st_point_on_surface st_polygonize st_precision
## [55] st_reverse st_sample st_segmentize
## [58] st_set_precision st_shift_longitude st_simplify
## [61] st_snap st_sym_difference st_transform
## [64] st_triangulate st_union st_voronoi
## [67] st_wrap_dateline st_write st_z_range
## [70] st_zm transform
## see '?methods' for accessing help and source code
st_crs(world)
## Coordinate Reference System:
## User input: WGS 84
## wkt:
## GEOGCRS["WGS 84",
## DATUM["World Geodetic System 1984",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## CS[ellipsoidal,2],
## AXIS["geodetic latitude (Lat)",north,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433]],
## AXIS["geodetic longitude (Lon)",east,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433]],
## USAGE[
## SCOPE["Horizontal component of 3D system."],
## AREA["World."],
## BBOX[-90,-180,90,180]],
## ID["EPSG",4326]]
st_crs(world)$IsGeographic
## [1] TRUE
st_crs(world)$units_gdal
## [1] "degree"
st_crs(world)$srid
## [1] "EPSG:4326"
st_crs(world)$proj4string
## [1] "+proj=longlat +datum=WGS84 +no_defs"
world_wgs84 = st_set_crs(world, "EPSG:4326") # set CRS to WGS84
st_crs(world_wgs84)$srid
## [1] "EPSG:4326"
Here we have a geographic coordinate system
#plot(world)
plot(world['continent'])
head(sf_proj_info(type = "proj"))
## name description
## 1 adams_hemi Adams Hemisphere in a Square
## 2 adams_ws1 Adams World in a Square I
## 3 adams_ws2 Adams World in a Square II
## 4 aea Albers Equal Area
## 5 aeqd Azimuthal Equidistant
## 6 affine Affine transformation
Here we have two projected coordinate systems
# Robinson Projections
world_robin <- st_transform(world, "+proj=robin")
plot(world_robin['lifeExp'])
# Mollweide Projections
world_moll <- st_transform(world, "+proj=moll")
plot(world_moll['gdpPercap'])
Here we can the the see the distortions in each of the systems for the USA
plot(world_robin['geom'][world_robin$name_long == 'United States', ])
plot(world_moll['geom'][world_moll$name_long == 'United States', ])
plot(world['geom'][world$name_long == 'United States', ])
We recommend to check out https://epsg.io for details about coordinate systems
Here we compare the US National Atlas Equal Area projection with WGS 84 coordinate systems
# the package data is stored in NAD83
st_crs(us_states)$srid
## [1] "EPSG:4269"
st_crs(us_states)
## Coordinate Reference System:
## User input: EPSG:4269
## wkt:
## GEOGCS["NAD83",
## DATUM["North_American_Datum_1983",
## SPHEROID["GRS 1980",6378137,298.257222101,
## AUTHORITY["EPSG","7019"]],
## TOWGS84[0,0,0,0,0,0,0],
## AUTHORITY["EPSG","6269"]],
## PRIMEM["Greenwich",0,
## AUTHORITY["EPSG","8901"]],
## UNIT["degree",0.0174532925199433,
## AUTHORITY["EPSG","9122"]],
## AUTHORITY["EPSG","4269"]]
# EPSG:2163 - US National Atlas Equal Area
p1 <- ggplot(data = us_states) +
geom_sf() +
coord_sf(crs = st_crs(2163)) + theme_bw()
# EPSG:4326 - WGS 84
p2 <- ggplot(data = us_states) +
geom_sf() +
coord_sf(crs = st_crs(4326)) + theme_bw()
library(patchwork)
p1 + p2
Test EPSG:2163 on world map
ggplot(data = world) +
geom_sf() +
coord_sf(crs = st_crs(2163)) + theme_bw()
Here we are combining two plots of Australia and the USA
us_sf <- world['geom'][world$name_long == 'United States', ]
us_sf <- st_transform(us_sf, "EPSG:4326")
aus_sf <- world['geom'][world$name_long == 'Australia', ]
aus_sf <- st_transform(aus_sf, "EPSG:4326")
us <- ggplot() +
geom_sf(data = us_sf)
aus <- ggplot() +
geom_sf(data = aus_sf)
us + aus
st_crs(us_sf)$srid
## [1] "EPSG:4326"
st_crs(aus_sf)$srid
## [1] "EPSG:4326"
Now we are plotting Australia and the USA in a single plot.
us_sf <- st_transform(us_sf, 'EPSG:4326')
aus_sf <- st_transform(aus_sf, 'EPSG:4326')
st_crs(us_sf)$units_gdal
## [1] "degree"
ggplot()+
geom_sf(data = us_sf)+
geom_sf(data = aus_sf)+
coord_sf(xlim = c(-130, 150), ylim = c(-50, 50), expand = FALSE) +
theme_void()
We can move Australia next to the USA. Keep in mind that we are still using WGS80 / EPSG:4326. That means we have distortions in our plot.
# copy the object
map_au_moved <- aus_sf
# transform (move) the object by -180 degrees + 60 degrees (north/west)
st_geometry(map_au_moved) <- st_geometry(map_au_moved) + c(-180, 60)
# copy the coordinate reference systems
st_crs(map_au_moved) <- st_crs(aus_sf)
ggplot()+
geom_sf(data = us_sf)+
geom_sf(data = map_au_moved)+
coord_sf(xlim = c(-130, -20), ylim = c(15, 50), expand = FALSE) +
theme_void()
## Warning in st_cast.GEOMETRYCOLLECTION(X[[i]], ...): only first part of
## geometrycollection is retained
To minimize distortion we can use local datums. Then we use these projections in coordinate system that reflect the same scale
us_sf <- st_transform(us_sf, 'EPSG:2163') # EPSG:2163 - US National Atlas Equal Area
aus_sf <- st_transform(aus_sf, 'EPSG:3112') # EPSG:3112 - GDA94 / Geoscience Australia Lambert
st_crs(us_sf)$units_gdal
## [1] "metre"
st_crs(aus_sf)$units_gdal
## [1] "metre"
ggplot()+
geom_sf(data = us_sf)+
# coordinate system with 6.000 km width and 4.000 km height
coord_sf(xlim = c(-3000000,3000000), ylim = c(-3000000,1000000)) +
theme_bw()
ggplot()+
geom_sf(data = aus_sf)+
# coordinate system with 6.000 km width and 4.000 km height
coord_sf(xlim = c(-3000000,3000000), ylim = c(-5000000,-1000000)) +
theme_bw()
Displaying a map of Europe by specifying the longitudes and latitudes with coord_sf() Notice: coord_sf() is applied on the whole world data.
ggplot() + geom_sf(data = world) +
coord_sf(xlim = c(-25, 50), ylim = c(30, 73), expand = FALSE) + theme_bw()
Here we filter our World dataset with ‘Europe’ The European overseas territories and Russia are not exactly what we wanted.
europe <-world['geom'][world$continent == 'Europe', ]
ggplot() + geom_sf(data = europe) + theme_bw()
Now we use the cropping function
st_crop() to cut our dataset.
europe_cropped <- st_crop(europe, xmin = -25, xmax = 50,
ymin = 30, ymax = 73)
ggplot() + geom_sf(data = europe_cropped) + theme_bw()
To locate the Hertie School on our map we first have to retrieve the coordinates. Then we can create a point geometry and plot it on top of our map.
coord_hertie <- c(13.389225005174465, 52.512803179279636)
hertie_school = st_sfc(st_point(coord_hertie), crs = 4326)
ggplot() + geom_sf(data = europe_cropped) +
geom_sf(data = hertie_school, color = 'red') +
theme_bw()
With st_buffer we can increase our geometry
hertie_buff = st_buffer(hertie_school, dist = 50 * 1000)
ggplot() + geom_sf(data = europe_cropped) +
geom_sf(data = hertie_buff, fill = 'red', alpha = 0.5, color = 'red') +
theme_bw()
The st_buffer also allows us to draw wider circles around our buffer. Since our CRS is still WGS80 we have distortions in Europe. Therefore, the circles are not perfect.
buff_500km = st_buffer(hertie_school, dist = 500 * 1000)
buff_1000km = st_buffer(hertie_school, dist = 1000 * 1000)
buff_1500km = st_buffer(hertie_school, dist = 1500 * 1000)
ggplot() + geom_sf(data = europe_cropped) +
geom_sf(data = hertie_buff, fill = 'red', alpha = 0.5, color = 'red') +
geom_sf(data = buff_500km, fill = 'grey', alpha = 0.1) +
geom_sf(data = buff_1000km, fill = 'grey', alpha = 0.1, ) +
geom_sf(data = buff_1500km, fill = 'grey', alpha = 0.1, ) +
theme_bw()
st_crs(europe_cropped)
## Coordinate Reference System:
## User input: WGS 84
## wkt:
## GEOGCRS["WGS 84",
## DATUM["World Geodetic System 1984",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## CS[ellipsoidal,2],
## AXIS["geodetic latitude (Lat)",north,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433]],
## AXIS["geodetic longitude (Lon)",east,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433]],
## USAGE[
## SCOPE["Horizontal component of 3D system."],
## AREA["World."],
## BBOX[-90,-180,90,180]],
## ID["EPSG",4326]]
To archive better circles we have to use a datum that is more suitable for Europe. For examples EPSG:5643 is recommended for central Europe.
new_europe <- st_transform(europe_cropped, 'EPSG:5643') # set CRS
new_buff <- st_transform(hertie_buff, 'EPSG:5643') # set CRS
buff_500km <- st_transform(buff_500km, 'EPSG:5643') # set CRS
buff_1000km <- st_transform(buff_1000km, 'EPSG:5643') # set CRS
buff_1500km <- st_transform(buff_1500km, 'EPSG:5643') # set CRS
ggplot() + geom_sf(data = new_europe) +
geom_sf(data = new_buff, fill = 'red', alpha = 0.5, color = 'red') +
geom_sf(data = buff_500km, fill = 'grey', alpha = 0.1) +
geom_sf(data = buff_1000km, fill = 'grey', alpha = 0.1, ) +
geom_sf(data = buff_1500km, fill = 'grey', alpha = 0.1, ) +
theme_bw()
st_crs(new_europe)$units_gdal
## [1] "metre"
st_crs(new_buff)$units_gdal
## [1] "metre"
st_crs(new_europe)$IsGeographic
## [1] FALSE
st_crs(new_buff)$IsGeographic
## [1] FALSE
Excercise 1
Try to retrive the spatial data for Canada from World country polygons (from the spData package)
Excercise 2
Try to find a suitable projection or datum for Canada on https://epsg.io and use it for plotting Canada in ggplot2.
Excercise 3
Try to find the coordinates for Canadas capital (Ottawa) and create a point geometry. Plot the point to that map you presivouly created.
Hopefully you learned something new! Feel free to contact us in case you have further questions and enjoy the rest of our interesting workshop day!